clc
clear all
close all 
addpath('functions');
rng(1234);

%% === Data === %% 
data.path           = '../B_temp/svar_data_final.csv'; 
dateformat          = 'yyyy''m''MM'; 
data.sample_y       = [datetime('1981-9-1'), datetime('2017-09-01')]; 
data.sample_z       = [datetime('1995-5-1'), datetime('2017-09-01')];  
data.y_names        = {'gs1', 'logcpi', 'logip', 'ebp','logTFP_interp','lmarkup_PF_w4_var_interp'};  
data.proxy_names    = 'ffr3shockTightWeighted';  

data.table = readtable(data.path); 

data.sample_y_bool = datetime(data.table.date) >= data.sample_y(1) & datetime(data.table.date) <= data.sample_y(2);
data.sample_z_bool = datetime(data.table.date) >= data.sample_z(1) & datetime(data.table.date) <= data.sample_z(2); 

data.time_y     = data.table.date(data.sample_y_bool); 
data.Y          = data.table{data.sample_y_bool, data.y_names}; 
data.time_z     = data.table.date(data.sample_z_bool); 
data.Z          = data.table{data.sample_z_bool, data.proxy_names}; 

VAR.data.Y      = data.Y;
VAR.data.time_y = data.time_y;
VAR.data.Z      = data.Z;
VAR.data.time_z = data.time_z; 

%% === Settings === %

settings.VAR_laglength  = 3;  
settings.IRF_hor = 49; 
settings.scale = 0.3;

VAR.p = settings.VAR_laglength;
VAR.IRF_hor = settings.IRF_hor;

col1 = [0 0 0];

%% === Estimation === %%
VAR = proxy_svar(VAR);    
VARbs = proxy_svar_bs(VAR, 10000, 68);    

%% === Plots === %%

for i=1:length(data.y_names) 
    h = figure;
    set(h,'DefaultAxesFontSize',20); 
    hold on
    hh=plot([0:length(VARbs.irs(:,i))-1],VARbs.irs(:,i)*0,'k-');
    plot([0:length(VARbs.irs(:,i))-1],VARbs.irs(:,i)/VARbs.irs(1,1)*settings.scale,'Color',col1,'Linewidth',6) 
    fan([0:length(VARbs.irs(:,i))-1],VARbs.irsH(:,i)/VARbs.irs(1,1)*settings.scale,VARbs.irsL(:,i)/VARbs.irs(1,1)*settings.scale,col1,'-',0.15,0);
        
    if ~strcmp(data.y_names{i},'lmarkup_PF_w4_var_interp') 
        ylabel('\%','Interpreter','Latex')
        hYLabel = get(gca,'YLabel');
        set(hYLabel,'rotation',0,'VerticalAlignment','middle')
        set(hYLabel, 'Units', 'Normalized', 'Position', [-0.15, 0.5, 0]);
    else
        ax = ancestor(hh, 'axes');
        ax.YAxis.Exponent = 0;
    end
    xlim([0 length(VARbs.irs(:,i))-1]) 
    xticks([0:12:48]) 
    xlabel('Months since shock','Interpreter','Latex')
    grid on
    box on
    hold off 
    fn_print(h,['''' '../C_results/VAR_' data.y_names{i} '.pdf''']); 
end